************************************************************************
* 				ADDITIONAL ROBUSTNESS AND CHECKS				       *
************************************************************************

use "data/iraqmerged.dta", clear


**********ISLAM ATTITUDES CHECKS*********

//Islam social desirability bias
local islvars = "isl1 isl2 isl4 isl5 isl6 isl7 isl8 isl9 isl10 isl11 isl12 isl13 isl14 isl15"
foreach var of local islvars {
	recode `var'(.=1)(1 2 3 4=0), into(`var'_mis) 
	}
	
local islmisvars = "isl1_mis isl2_mis isl4_mis isl5_mis isl6_mis isl7_mis isl8_mis isl9_mis isl10_mis isl11_mis isl12_mis isl13_mis isl14_mis isl15_mis"
		
	foreach var of local islmisvars {
	logit `var' treat i.regid if kurd==0 &shia==0 
	}
*Sunnis significantly less likely to be missing on isl10_mis and isl5_mis and isl4_mis, significantly more likely for isl7_mis

local islmisvars = "isl1_mis isl2_mis isl4_mis isl5_mis isl6_mis isl7_mis isl8_mis isl9_mis isl10_mis isl11_mis isl12_mis isl13_mis isl14_mis isl15_mis"	
	foreach var of local islmisvars {
	logit `var' treat i.regid   if kurd==0 &shia==1
	}
*Shi'is significantly more likely to be missing on isl15_mis isl13_mis isl6_mis


//Check rejection of political Islam
regress relmean treat##shia
mixed relmean treat##shia c.age i.gender i.coledu c.pctshia c.sqsect c.clanperthou if kurd==0 || regid: || disid:

//Check for any sectarian dimension of Arab identity
melogit arabident c.age i.gender i.coledu c.pctshia c.sqsect c.clanperthou c.relmean if kurd==0 & shia==0 & treat==1 || regid: || disid:

//Is effect stronger for people who consume more news?
melogit natident i.treat##i.shia i.treat##c.newsint c.age i.gender i.coledu c.clanperthou c.sqsect c.pctshia if kurd==0 [pw=_webal] || regid: || disid:

//Are news readers more worried about war?
mixed secint i.treat##i.shia i.treat##c.newsint c.age i.gender i.coledu c.clanperthou c.sqsect c.pctshia if kurd==0 || regid: || disid:
mixed civwarint i.treat##i.shia i.treat##c.newsint c.age i.gender i.coledu c.clanperthou c.sqsect c.pctshia if kurd==0 || regid: || disid:
mixed warint i.treat##i.shia i.treat##c.newsint c.age i.gender i.coledu c.clanperthou c.sqsect c.pctshia if kurd==0 || regid: || disid:

**********TIME TRENDS*********

//Generate time running variable where 0 corresponds to the day of Mosul's fall
gen time_zero =datevar-19884

local idents = "natident arabident islident"
foreach outvar of local idents{
	preserve
	keep if kurd==0 &shia==1
	bysort time_zero: egen avg`outvar' = mean(`outvar')
	graph twoway (scatter avg`outvar' time_zero if time_zero < 0,  ysc(r(0 1))  ylabel(0(.5)1) xlabel(-12(4)0) ) ///
				(lpoly avg`outvar' time_zero if time_zero < 0) /// 
				(lowess avg`outvar' time_zero if time_zero < 0) ///
				(lfit avg`outvar' time_zero if time_zero < 0, lpattern(solid)), /// 
				xtitle("Days to event") ylab(,nogrid) scheme(s1mono) ///
				legend(order(2 "Kernel weighted local polynomial smoothing" 3 "Locally weighted regression smoothing" 4 "Linear Fit") size(vsmall) rows(2) pos(6)) ///
				ysize(3) name(ttshi`outvar', replace)
restore
}

foreach outvar of local idents{
	preserve
	keep if kurd==0 &shia==0
	bysort time_zero: egen avg`outvar' = mean(`outvar')
	graph twoway (scatter avg`outvar' time_zero if time_zero < 0,  ysc(r(0 1))  ylabel(0(.5)1) xlabel(-12(4)0) ) ///
				(lpoly avg`outvar' time_zero if time_zero < 0) /// 
				(lowess avg`outvar' time_zero if time_zero < 0) ///
				(lfit avg`outvar' time_zero if time_zero < 0, lpattern(solid)), /// 
				xtitle("Days to event") ylab(,nogrid) scheme(s1mono) ///
				legend(order(2 "Kernel weighted local polynomial smoothing" 3 "Locally weighted regression smoothing" 4 "Linear Fit") size(vsmall) rows(2) pos(6)) ///
				ysize(3) name(ttsun`outvar', replace)
restore
}


grc1leg ttshinatident ttsunnatident ttshiarabident ttsunarabident ttshiislident ttsunislident, rows(3)
gr_edit plotregion1.graph1.yaxis1.title.text.Arrpush Shi'i Nat ident.
gr_edit plotregion1.graph2.yaxis1.title.text.Arrpush Sunni Nat ident.
gr_edit plotregion1.graph3.yaxis1.title.text.Arrpush Shi'i Arab ident.
gr_edit plotregion1.graph4.yaxis1.title.text.Arrpush Sunni Arab ident.
gr_edit plotregion1.graph5.yaxis1.title.text.Arrpush Shi'i Islamic ident.
gr_edit plotregion1.graph6.yaxis1.title.text.Arrpush Sunni Islamic ident.

graph export "plots/ttidents.png", replace


//Generate combined security var.
gen secallint = secint + secregint + secnatint

local secvars = "secallint secint secregint secnatint"


foreach outvar of local secvars {
	bysort time_zero: egen avg`outvar' = mean(`outvar') if kurd==0
	graph twoway (scatter avg`outvar' time_zero if time_zero < 0,  ysc(r(0 12))  ylabel(0(5)15) xlabel(-12(4)0) ) ///
				(lpoly avg`outvar' time_zero if time_zero < 0) /// 
				(lowess avg`outvar' time_zero if time_zero < 0) ///
				(lfit avg`outvar' time_zero if time_zero < 0, lpattern(solid)), ///
				xtitle("Days to event") ylab(,nogrid) scheme(s1mono) ///
				legend(order(2 "Kernel weighted local polynomial smoothing" 3 "Locally weighted regression smoothing" 4 "Linear Fit") size(vsmall) rows(2) pos(6)) ///
				ysize(3) name(sec_`outvar', replace)
}

grc1leg sec_secallint sec_secint sec_secregint sec_secnatint
gr_edit plotregion1.graph1.yaxis1.title.text.Arrpush Combined security rating
gr_edit plotregion1.graph2.yaxis1.title.text.Arrpush Local security rating
gr_edit plotregion1.graph3.yaxis1.title.text.Arrpush Regional security rating
gr_edit plotregion1.graph4.yaxis1.title.text.Arrpush National security rating

graph export "plots/ttsecints.png", replace

local warvars = "warint civwarint"

foreach outvar of local warvars {
	bysort time_zero: egen avg`outvar' = mean(`outvar') if kurd==0
	graph twoway (scatter avg`outvar' time_zero if time_zero < 0,  ysc(r(0 4))  ylabel(0(1)4) xlabel(-12(4)0)) ///
				(lpoly avg`outvar' time_zero if time_zero < 0) /// 
				(lowess avg`outvar' time_zero if time_zero < 0) ///
				(lfit avg`outvar' time_zero if time_zero < 0, lpattern(solid)), ///
				xtitle("Days to event") ylab(,nogrid) scheme(s1mono) ///
				legend(order(2 "Kernel weighted local polynomial smoothing" 3 "Locally weighted regression smoothing" 4 "Linear Fit") size(vsmall) rows(2) pos(6)) ///
				ysize(3) name(war_`outvar', replace) aspectratio(.3)
}

grc1leg war_warint  war_civwarint, rows(2)
gr_edit plotregion1.graph1.yaxis1.title.text.Arrpush War worry
gr_edit plotregion1.graph2.yaxis1.title.text.Arrpush Civil war worry

graph export "plots/ttwarints.png", replace


**********STAT POWER BANDWIDTH TESTS*********
drop time_zero
encode date1, gen(decdate)
//Pre-existing time trends
gen yearst=substr(date1,1,4)
gen monthst=substr(date1,5,2)
gen dayst=substr(date1,7,2)
destring yearst, gen(year)
destring monthst, gen(month)
destring dayst, gen(day)
gen edate = mdy(month, day, year)
gen timepoint = edate - 19870
gen time_zero = decdate - 9
replace time_zero = . if decdate == 9
replace time_zero = time_zero - 1 if time_zero > 0

gen treatment_naive=0
replace treatment_naive=1 if time_zero>=0
replace treatment_naive=. if time_zero==.


* Generating matrix that includes the power and number of units in the treatment and control groups for two differnt effect sizes and multiple bandwidths
preserve 
drop if kurd==1|shia==0
matrix results = J(20,5,.)
matrix colnames results =  bandwidth n_contr n_treat pow_thrsd pow_halfsd
local i 0 
forval d = 0/8{
	if `d' == 0{
		local d_neg = -1
	}
	if `d' != 0{
				local d_neg = (`d'+1) * -1 
	}
	gen treatment_band = . 
	replace treatment_band = 1 if time_zero >= 0 & time_zero <= `d' 
	replace treatment_band = 0 if time_zero < 0 & time_zero >= `d_neg'
	replace treatment_band = . if missing(age, gender, coledu)								 
	local ++ i 
	local bandwidth = `d' + 1 
	qui: sum treatment_band if treatment_band == 1 
	local n_treat = r(N)
	qui: sum treatment_band if treatment_band == 0 
	local n_contr = r(N)
	sum natident if treatment_naive == 0 
	local sd = r(sd) 
	local mean = r(mean)
	local half = `sd' / 2
	local thir = `sd' / 3
	local thirsd = `mean' + `thir'
	local halfsd = `mean' + `half'
	power twomeans `mean' `thirsd', sd(`sd') n1(`n_contr') n2(`n_treat') 
	local pow_thrsd = r(power)
	power twomeans `mean' `halfsd', sd(`sd') n1(`n_contr') n2(`n_treat')
	local pow_halfsd = r(power)	
	matrix results[`i', 1] = `bandwidth',`n_contr',`n_treat',`pow_thrsd',`pow_halfsd'
	drop treatment_band
} 

* Transforming matrix into variables
svmat results, names(col)

* Total height of histogram bars
gen tot_contr = n_treat + n_contr

* Generating Figure 
local plus = char(177)

graph twoway	(bar tot_contr bandwidth, fcolor(gs4) barw(0.85) ///
				yaxis(2) ylab(,nogrid) ytitle("Number of valid cases", axis(2)) ///
				yscale(axis(2) alt)) ///
				(rbar tot_contr n_treat bandwidth, fcolor(gs13) barw(0.85) yaxis(2)) ///
				(function y=.8, ra(0 10) lpattern(solid) lcolor(black) lwidth(vthin)) ///
				(line pow_thrsd bandwidth, lcolor(black) lpattern(solid) lcolor(black) yaxis(1) yscale(axis(1) alt)) ///
				(line pow_halfsd bandwidth, lpattern(dash) lcolor(black) yaxis(1) ///
				ytitle("Power") ylabel(0(.2)1, gmin gmax axis(1)) ///
				xtitle("Bandwidth (± days) around event") ///
				scheme(s1mono) ///
				legend(ti("Effect size", size(vsmall)) ///
				order(1 "Control group" 2 "Treatment group" 4 "1/3 standard deviation change" 5 "1/2 standard deviation change") size(vsmall) rows(2) pos(6)))
				
gr_edit style.editstyle boxstyle(shadestyle(color(white))) editcopy
gr_edit .plotregion1.plot2.style.editstyle area(linestyle(width(none))) editcopy
gr_edit .plotregion1.plot2.style.editstyle area(linestyle(width(none))) editcopy

graph export "plots/statpow1.png", replace
restore


preserve 

drop if kurd==1|shia==1
matrix results = J(20,5,.)
matrix colnames results =  bandwidth n_contr n_treat pow_thrsd pow_halfsd
local i 0 
forval d = 0/8{
	if `d' == 0{
		local d_neg = -1
	}
	if `d' != 0{
				local d_neg = (`d'+1) * -1 
	}
	gen treatment_band = . 
	replace treatment_band = 1 if time_zero >= 0 & time_zero <= `d' 
	replace treatment_band = 0 if time_zero < 0 & time_zero >= `d_neg'
	replace treatment_band = . if missing(age, gender, coledu)								 
	local ++ i 
	local bandwidth = `d' + 1 
	qui: sum treatment_band if treatment_band == 1 
	local n_treat = r(N)
	qui: sum treatment_band if treatment_band == 0 
	local n_contr = r(N)
	sum arabident if treatment_naive == 0 
	local sd = r(sd) 
	local mean = r(mean)
	local half = `sd' / 2
	local thir = `sd' / 3
	local thirsd = `mean' + `thir'
	local halfsd = `mean' + `half'
	power twomeans `mean' `thirsd', sd(`sd') n1(`n_contr') n2(`n_treat') 
	local pow_thrsd = r(power)
	power twomeans `mean' `halfsd', sd(`sd') n1(`n_contr') n2(`n_treat')
	local pow_halfsd = r(power)	
	matrix results[`i', 1] = `bandwidth',`n_contr',`n_treat',`pow_thrsd',`pow_halfsd'
	drop treatment_band
} 

* Transforming matrix into variables
svmat results, names(col)

* Total height of histogram bars
gen tot_contr = n_treat + n_contr

* Generating Figure 
local plus = char(177)

graph twoway	(bar tot_contr bandwidth, fcolor(gs4) barw(0.85) ///
				yaxis(2) ylab(,nogrid) ytitle("Number of valid cases", axis(2)) ///
				yscale(axis(2) alt)) ///
				(rbar tot_contr n_treat bandwidth, fcolor(gs13) barw(0.85) yaxis(2)) ///
				(function y=.8, ra(0 10) lpattern(solid) lcolor(black) lwidth(vthin)) ///
				(line pow_thrsd bandwidth, lcolor(black) lpattern(solid) lcolor(black) yaxis(1) yscale(axis(1) alt)) ///
				(line pow_halfsd bandwidth, lpattern(dash) lcolor(black) yaxis(1) ///
				ytitle("Power") ylabel(0(.2)1, gmin gmax axis(1)) ///
				xtitle("Bandwidth (± days) around event") ///
				scheme(tufte) ///
				legend(ti("Effect size", size(vsmall)) ///
				order(1 "Control group" 2 "Treatment group" 4 "1/3 standard deviation change" 5 "1/2 standard deviation change") size(vsmall) rows(2) pos(6)))
				
gr_edit style.editstyle boxstyle(shadestyle(color(white))) editcopy
gr_edit .plotregion1.plot2.style.editstyle area(linestyle(width(none))) editcopy
gr_edit .plotregion1.plot2.style.editstyle area(linestyle(width(none))) editcopy

graph export "plots/statpow2.png", replace
restore


*******T-TEST BALANCE PLOTS BEFORE AND AFTER******

use "data/iraqmerged.dta", clear

gen treatment_0 = treat
gen treatment_1 = .
replace treatment_1=0 if DATEINT<20140610
replace treatment_1=1 if DATEINT>20140610

* Reversing treatment variables so that in means comparisons higher numbers means higher (proportion) in treatment group
foreach var of varlist treatment_*{
	gen r_`var' = `var'
	recode r_`var' (1=0) (0=1)
}
		
		
*global indvarsecvstd3 r_age gender coledu wealthlm wealthmh wealthh unemp housewife retired inwork voted hajj rels rel r_relmean gentrust
*global disvarsecvstd3 r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist


			
local varlist age pctshia sqmosdist sqsect popdens clanperthou nmis intlsecs
foreach v of var `varlist' {
qui summarize `v'
gen r_`v' = `v'
replace r_`v' = r_`v'/ (r(sd) * 2)
}

			
foreach tr of varlist r_treatment_*{
	matrix mean = J(1,28,.)
	matrix colnames mean = r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist
	matrix CI = J(4,28,.)
	matrix colnames CI = r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist
	matrix rownames CI = ll95 ul95 ll90 ul90
	local i 0
	foreach var of varlist r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist {
		quietly: ttest `var', by(`tr') 
		local ++ i 
		local diff = r(mu_1) - r(mu_2) 
		matrix mean[1, `i'] = `diff' 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
		matrix CI[1, `i'] = `ll95' \ `ul95' \ `ll90' \ `ul90'
	}
matrix `tr'_m = mean
matrix `tr'_CI = CI
}


* Generating Figure from results stored in matrices 
local plus = char(177)

lab var r_age "Age"
lab var r_pctshia "% Shi'i'"
lab var r_sqmosdist "Dist. from Mosul"
lab var r_sqsect "Prehist. violence"
lab var r_clanperthou "Tribal ties"
lab var urban "Urban/rural"
lab var r_popdens "Pop. density"
lab var r_nmis "Missing responses"
lab var r_intlsecs "Completion time"

preserve
keep if shia==0
foreach tr of varlist r_treatment_*{
	matrix mean = J(1,28,.)
	matrix colnames mean = r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist
	matrix CI = J(4,28,.)
	matrix colnames CI = r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist
	matrix rownames CI = ll95 ul95 ll90 ul90
	local i 0
	foreach var of varlist r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist {
		quietly: ttest `var', by(`tr') 
		local ++ i 
		local diff = r(mu_1) - r(mu_2) 
		matrix mean[1, `i'] = `diff' 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
		matrix CI[1, `i'] = `ll95' \ `ul95' \ `ll90' \ `ul90'
	}
matrix `tr'_msun = mean
matrix `tr'_CIsun = CI
}
restore


preserve
keep if shia==1
foreach tr of varlist r_treatment_*{
	matrix mean = J(1,28,.)
	matrix colnames mean = r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist
	matrix CI = J(4,28,.)
	matrix colnames CI = r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist
	matrix rownames CI = ll95 ul95 ll90 ul90
	local i 0
	foreach var of varlist r_age gender coledu wealthl wealthlm wealthmh wealthh  voted hajj  r_nmis r_intlsecs r_pctshia r_sqsect r_clanperthou r_popdens urban r_sqmosdist {
		quietly: ttest `var', by(`tr') 
		local ++ i 
		local diff = r(mu_1) - r(mu_2) 
		matrix mean[1, `i'] = `diff' 
		local degrees = r(df_t)
		local critical_5 = invttail(`degrees', 0.025)
		local confvalue_5 = `critical_5' * r(se)
		local critical_10 = invttail(`degrees', 0.05)
		local confvalue_10 = `critical_10' * r(se)
		local ll95 = `diff' - `confvalue_5'
		local ul95 = `diff' + `confvalue_5'
		local ll90 = `diff' - `confvalue_10'
		local ul90 = `diff' + `confvalue_10'
		matrix CI[1, `i'] = `ll95' \ `ul95' \ `ll90' \ `ul90'
	}
matrix `tr'_mshi = mean
matrix `tr'_CIshi = CI
}
restore

coefplot	(matrix(r_treatment_0_m), xline(0, lwidth(vvthin))	///
			ci((r_treatment_0_CI[1] r_treatment_0_CI[2]))) ///
			|| (matrix(r_treatment_0_mshi), xline(0, lwidth(vvthin))	///
			ci((r_treatment_0_CIshi[1] r_treatment_0_CIshi[2]) )) ///
			|| (matrix(r_treatment_0_msun), xline(0, lwidth(vvthin))	///
			ci((r_treatment_0_CIsun[1] r_treatment_0_CIsun[2]) )) ///
			, byopts(row(1)) aspectratio(1) scheme(s1mono) ciopts(lwidth(thin)) msize(vsmall) ///
			xlab(,nogrid labsize(vsmall)) ylab(, labsize(vsmall)) ytick(11.5, notick glstyle(refline) glpattern(dash) glcolor(maroon) glwidth(.3)) ///
			bylabels("Full sample" "Shi'i" "Sunni") subtitle(, size(small))

graph export "plots/ttestsall.pdf", replace

**********PLACEBO TESTS*********

use "data/iraqmerged.dta", clear
encode date1, gen(decdate)
//Pre-existing time trends
gen yearst=substr(date1,1,4)
gen monthst=substr(date1,5,2)
gen dayst=substr(date1,7,2)
destring yearst, gen(year)
destring monthst, gen(month)
destring dayst, gen(day)
gen edate = mdy(month, day, year)
gen timepoint = edate - 19870
gen time_zero = decdate - 9

gen time_zerop = time_zero
replace time_zerop = . if decdate == 9
replace time_zerop = time_zerop - 1 if time_zerop > 0

gen treatment_naive=0
replace treatment_naive=1 if time_zerop>=0
replace treatment_naive=. if time_zerop==.

/*
Generating placebo treatment variable at the empirical median of the control group (June 6) for both sects.
Note June 6 is coded as 1 for placebo treatment. Same results obtain when coded as . (in line with suggestion
of Munoz to code day of event as missing).
*/

gen placebo = . 
sum decdate if treatment_naive == 0 , det
replace placebo = 0 if decdate <5
replace placebo = 1 if decdate >=5 & decdate <=8
tab placebo
tab decdate placebo, m
gen time_interview = decdate
label var placebo "Placebo treatment"
label var time_interview "Fieldwork day"

global placsect i.placebo##i.shia
global indvars c.age i.gender i.coledu
global disvars c.pctshia c.sqsect c.clanperthou

melogit natident $placsect $indvars $disvars if kurd==0 || regid: || disid: 
est store pmn
melogit arabident $placsect $indvars $disvars if kurd==0 || regid: || disid: 
est store pma

local modlist = "pmn pma"

foreach model of local modlist {

//make separate predictions for shia and sunni at treat=0/1
est resto `model'
margins r.placebo, over(r.shia) noestimcheck post
est sto mc`model'
}

coefplot (mcpmn), bylabel({bf:Nat. ident.}) || ///
(mcpma), bylabel({bf:Arab. ident.}) ||, ///
vertical yline(0, lcolor(black) lwidth(vvthin)) title(,size(tiny)) ytitle("Contrasts, Shi'i : Sunni", size(vsmall)) title("") ///
		levels(95) ylab(,nogrid labsize(vsmall)) xlab("",nogrid labsize(vsmall)) ciopts(lwidth(thin)) ///
		scheme(s1mono) ///
		subtitle(, size(small)) ///
		mlabposition(2) mlabsize(vsmall) ///
		mlabel format(%12.2f) mlabposition(3.5) mlabgap(*2) mlabsize(vsmall) ///
		legend(rows(1) size(small)) byopts(row(1)) title("Model specification", size(3.5)) aspectratio(.3)
gr_edit plotregion1.title[1].text = {}
gr_edit plotregion1.title[2].text = {}

graph export "plots/MLmcspt.pdf", replace


//Placebo using Jordan ArabTrans_data
use "daya/raw//Arab Transformations/Arab Transformations Data Set STATA Version PA20170503.dta", clear
keep if COUNTRY==4

decode DATEINT, gen(date1)
gen datevar=date(date1,"YMD")
*br DATEINT date1 datevar
gen double longdate = date(date1,"YMD")
format longdate %td
encode date1, gen(decdate)
gen treat = 0
replace treat=1 if DATEINT>=20140610
recode V28(1=1 "Nation")(2 3 4 5=0 "Other")(6 98=.), into(natident)
recode V28(2=1 "Arab")(1 3 4 5=0 "Other")(6 98=.), into(arabident)
gen age=2014-V76
recode V77(1=0 "Male")(2=1 "Female"), into(gender)
gen coledu=0
replace coledu=1 if V79==6|V79==7
rename DISTRICT disid
rename REGION regid

melogit natident i.treat age gender coledu || regid: || disid:
est store jordpn

melogit arabident i.treat age gender coledu || regid: || disid:
est store jordpa


coefplot (jordpn), bylabel({bf:Nat. ident.}) || ///
jordpa, bylabel({bf:Arab. ident.}) ||, drop(_cons age gender coledu) vertical yline(0, lcolor(black) lwidth(vvthin)) title(,size(tiny)) ytitle("Coefficient on June 10 dummy", size(vsmall)) ///
		levels(95) ylab(,nogrid labsize(vsmall)) xlab("",nogrid labsize(vsmall)) ciopts(lwidth(thin)) ///
		scheme(s1mono) ///
		subtitle(, size(small)) ///
		mlabposition(2) mlabsize(vsmall) ///
		mlabel format(%12.2f) mlabposition(3) mlabgap(*2) mlabsize(vsmall) ///
		legend(rows(1) size(small) lcolor(none)) byopts(row(1)) name(MLptjor, replace) title("", size(5)) aspectratio(.3)

graph export "plots/MLptjor.pdf", replace


//Descriptives

global treatcovintsall i.treat#c.age i.treat#i.gender i.treat#i.coledu i.treat#c.pctshia i.treat#c.sqsect
global indvarsall c.age i.gender i.coledu i.wealthcatr i.empstat i.voted i.hajjbin
global disvarsall c.pctshia c.sqsect c.clanperthou c.popdens i.urban c.sqmosdist
global rqualvars c.nmis c.intlsecs

eststo clear
estpost sum natident arabident islident age gender coledu wealthcatr voted hajjbin pctshia sqsect clanperthou popdens urban sqmosdist nmis intlsecs if kurd==0
est store a
 
esttab a using latex_iraq/regtabs_v2/destab.tex, replace ///
collabels(\multicolumn{1}{c}{{Mean}} \multicolumn{1}{c}{{Std.Dev.}} \multicolumn{1}{l}{{Obs}}) ///
cells("mean(fmt(2)) sd(fmt(2)) count(fmt(0))") label nonumber f noobs alignment(S) booktabs


//Interflex

set scheme s1mono
//interflex natident treat clanperthou age gender coledu pctshia sqsect if kurd==0 &shia==0, ylab(Nat ident.) dlab(Treat) xlab(Tribal ties) 
//gr_edit plotregion1.plot4.bar_size = 3
//gr_edit plotregion1.plot3.bar_size = 3
interflex natident treat clanperthou age gender coledu pctshia sqsect if kurd==0 &shia==0, ylab(Nat ident.) dlab(Treat) xlab(Tribal ties) xd(density)

//Comparison with MEVS
use "data/raw/FinalDataCollectionDataset__CrossNationalStudyofValuesandValuesChange.dta", clear
keep if country==368
recode v58(1=1 "National identity")(2 3 4=0 "Other"), into(natident)
****high tribal ties crosstabs

//Salah ad-Din
tab natident if shia==0 & kurd==0 & regid == 40006

//Diyala
tab natident if shia==0 & kurd==0 & regid == 40007
